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Abstract 

We consider two-layers of immiscible liquids confined between an upper and a lower rigid plate. 
The dynamics of the free liquid-liquid interface is described for arbitrary amplitudes by a single 
evolution equation derived from the basic hydrodynamic equations using long-wave approximation. 
After giving the evolution equation in a general way, we focus on interface instabilities driven by 
gravity, thermocapillary and electrostatic fields. First, we study the linear stability discussing 
especially the conditions for destabilizing the system by heating from above or below. Second, we 
use a variational formulation of the evolution equation based on an energy functional to predict 
metastable states and the long-time pattern morphology (holes, drops or maze structures). Finally, 
fully nonlinear three-dimensional numerical integrations are performed to study the short- and long- 
time evolution of the evolving patterns. Different coarsening modes are discussed and long-time 
scaling exponents are extracted. 



2 



I. INTRODUCTION 



The study of stability properties and pattern formation in thin films is still an important 
and not fully explored challenge in fluid dynamics. Reorganization processes of such films 
have remarkable applications in chemical engineering, biological systems, or semiconductor 
industry. Besides industrial aspects the computational advantage of studying thin films is 
obvious since one equation for the surface is often sufficient to capture the basic dynamics. 
Due to the increase of computer power, pattern formation in systems far from equilibrium can 
be investigated in more detail. This leads to consideration of more and more complex systems 
which may show a rich variety of bifurcations and allows for a more realistic modelling of 
fluid phenomena. 

Here, we consider thin two-layer films bounded by an upper and a lower rigid plate. Using 
lubrication approximation a general long-wave interface evolution equation is derived that is 
valid for arbitrary amplitudes of interface deflections. Pattern formation under the influence 
of gravity, thermocapillarity and electrostatic forces is analysed. 

Lubrication or long- wave approximation is used for more than 100 years to describe the 
evolution of thin filmsi. In one-layer systems with a free interface the dynamics of the sur- 
rounding gas is normally neglected and solely the liquid determines the interface evolution. 
A simplified equation for the evolution of the profile of the surface can be derived from 
the basic hydrodynamic equations because the velocity is enslaved to the thickness profile. 
Several mechanisms are known to destabilize an initially flat surface. A survey of long-wave 
instabilities in one-layer systems is given by Oron et al.^. A prominent example is the 
destabilization and subsequent evolution of a non-flat surface profile due to Marangoni flow 
caused by heating from below. It was first studied by SCRIVEN and Sterling^ and classified 
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by Goussis and Kelly^ as the S-mode instability. The second mode, called P-mode, is a 
short-wave instability without surface deflection and will not concern us here. However, see 
GOLOVIN ET AL.^ for a study of the interaction between short- and long-wave mode. The 
linear and nonlinear behavior for an unstable thin liquid layer heated from below is studied 
by BuRELBACH ET AL.-. Deissler and Oron'' show the stabilizing effect of thermocap- 
illarity for a thin film at the underside of a cooled horizontal plate which is gravitationally 
or Rayleigh- Taylor (RT) unstable. The normally used linear dependence of surface ten- 
sion on temperature (linear Marangoni effect) is replaced by Oron and Rosenau^ by a 
quadratic one, thereby inhibiting true film rupture. Three-dimensional simulations of the 
linear Marangoni effect are done by Oron^ and for a wetting film by Bestehorn et AL.— . 
The former work concentrates on the evolution towards film rupture whereas the latter sys- 
tem allows to explain the preference of drops or holes in dependence of the film thickness. 
It also gives scaling exponents for the long-time coarsening. In a recent work Thiele and 
KNOBLOCHii show that the rich bifurcation structure for drop solutions on a horizontal 
substrate is destroyed even by a rather small inclination of the substrate. 
Mathematically related to thin heated films are ultrathin free surface films on horizontal 
substrates as first studied in long-wave approximation by Ruckenstein and jAlls^i^. These 
films may be unstable and dewet due to effective molecular interactions that are incorporated 
into the governing equations in form of an additional pressure term. This so-called disjoining 
pressure was introduced by Derjaguin et al. (for an overview see^^). In the simplest 
case it results from the apolar London-van der Waals dispersion forces^^. Open questions 
regarding dewetting in one-layer systems are summarized in Ref.— . Here, we will use a 
stabilizing van der Waals interaction to avoid 'true' film rupture for a heated filn>i^. 
The evolution of unstable thin films has a general characteristics that is found as well for thin 
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heated films as also for ultrathin dewetting films. One distinguishes between a short-time 
and a long-time behavior. First, the fiat film evolves into a large amplitude pattern whose 
typical length scale can normally be determined by linear considerations. Often, this stage 
is called initial film rupture although the film may not rupture totally, but an ultrathin film 
remains at the 'dry' parts. The long-time behavior is characterized by an ongoing coarsening 
towards patterns of longer and longer horizontal spatial scales. 

Evidently, long-wave approximation is not only applicable for single liquid layers on a solid 
substrate. The approach can be naturally extended towards systems characterized by more 
than one layer. Taking, for example, two layers of immiscible liquids on a solid horizontal 
substrate in a gas atmosphere yields a pair of coupled evolution equations for the liquid-liquid 
and the liquid-gas interfaces. Such a system in presence of a surfactant and an evaporating 
upper liquid is considered by Danov et AhMi^^^. Different pathways of dewetting induced 
by long-range van der Waals forces are investigated by Pototsky et al.^^. However, if 
such a two-layer system is bounded below and above by rigid plates its behavior can be 
described by a single evolution equation for the liquid-liquid interface. This kind of system 
is treated in the present paper. 

Although a general evolution equation was to our knowledge not yet given in the literature 
there exist a number of analyses for related systems. Yiantsios and HlGGlNS^^ investigate 
the Ray leigh- Taylor instability regarding an upper layer of infinite thickness. They use 
lubrication approximation for the lower liquid layer but not for the infinitely extended upper 
one. They find that to leading order the dynamics of the upper layer can be neglected if the 
viscosities of both liquids are of the same order of magnitude. In this way, they obtain an 
effective one-layer interface evolution equation. 

Marangoni effects in two superposed fluid layers are experimentally studied by VanHook 
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ET AL.^°. They investigate long-wave as well as short-wave thermocapillary instabilities. 
However, their theoretical analysis neglects velocities in the upper layer and uses a 'two-layer 
Biot number' to take into account the thermal properties as well as the thermal field in the 
upper layer. This generalization of the Biot number used in2. also leads to an effective one- 
layer equation. A similar theory is used by BuRGESS ET AL.^^ to explain the stabilization of 
a Rayleigh- Taylor (RT) unstable oil-air system by heating form below. Linear investigations 
of short- and long- wave Marangoni instabilities in two superposed liquid films are presented 
by Smitb2^. Furthermore, SiMANOVSKll and NEPOMNYASHCHYi^^ consider a two-layer 
system with thermocapillary effects. They derive a weakly nonlinear interface equation 
in long-wave approximation taking into account the dynamics in both layers. Our linear 
results for the thermocapillary case can be directly compared to theirs. They show that 
the occurrence of thermocapillary instabilities is not only determined by the direction of 
the temperature gradient but also by the ratios of the viscosities and the layer thicknesses. 
In particular, they find that contrary to an one-layer system heating from above can act 
destabilizing. Moreover, Tilley et. Ab^^ investigate two superposed fiuids in an inclined 
channel with gravity and Marangoni effects. Their weakly nonlinear analysis reveals a 
modified Kuramoto-Sivashinsky equation with broken refiectional symmetry. 
Two superposed dielectric fiuids between two parallel plates are an appropriate system to 
investigate pattern formation through electrohydrodynamic instabilities since a vertically 
applied electric field causes normal and tangential interface forces which depend strongly 
on the dielectric fiuid properties. Majumdar and O'NeilIj^S propose an experimental 
method to quantify surface tension via the measured critical voltage for the onset of such 
an instability. Mohamed et al.^^ investigate the destabilization of the interface using 
an Orr-Sommerfeld equation. Our linear results can be compared to theirs for quadratic 
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velocity profiles since long-wave approximation allows for quadratic velocity profiles only. A 
detailed analysis of different electrical fluid properties like the creation of free charges at the 
interface, or the characterization of conducting and insulating fluids is given by Melcher 
and Smithes,. Investigations of a free thin liquid jet with an axial applied electric field are 
done by Savettaseranee et al.^^. They show that the electric field stabilizes the film 
and can avoid rupture induced by attractive van der Waals forces. Experiments of LiN ET 
AL.^Si^ using two layers of polymeric liquid exposed to a vertical electric field suggest that 
linear considerations do well capture the length scale found even in the long-time evolution. 
Our nonlinear calculations confirm the validity of the linear theories. 

The present work focuses on two-layer films bounded by an upper and a lower rigid plate 
as sketched in Fig.Q In Section |n] we derive from the basic hydrodynamic equations the 
interface evolution equation in lubrication or long-wave approximation for general layer 
properties. Keeping the normal and tangential stresses at the interface in a general form, 
the resulting equation can be applied to the study of various body and interface forces. In 
the remaining paper we focus on gravity, thermo capillarity and electrostatic forces. This 
allows for the formulation of the problem in variational form using a Lyapunov functional. 
Since the free energy density is a function of the interface h only, the long-time evolution 
can be predicted, i.e. whether holes, drops or maze structures are expected for t oo. 
Further on we discuss the occurrence of metastable states. In Section IIIH we perform 
linear and nonlinear analyses of the derived equation. First, we show that gravitation and 
thermocapillarity may act stabilizing as well as destabilizing depending on material and 
system parameters. Furthermore, we integrate the fully nonlinear evolution equation in 
three dimensions numerically and study the short-time as well as the long-time evolution. 
For the latter, different coarsening modes and the long-time scaling are discussed. We 
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summarize our results in Section |TV| and point out possible applications, in particular the 
influence of electrostatic fields on dewetting. In the Appendix we discuss shortly the subtle 
occurrence of an additional mean flow field if the system is extended from 2D to 3D. 

II. GOVERNING EQUATIONS 

We consider a two-layer system bounded by a rigid upper and lower plate with a system 
thickness d and a fiat film interface with height (Fig.Q). Instabilities lead to a time and 
space dependent interface profile h{x,y,t). 

A. Evolution equation 

The required material parameters for incompressible fluids are the densities pi and the 
dynamic viscosities /ij. The subscripts i = 1 and i = 2 denote liquid 1 (lower layer) and 
liquid 2 (upper layer), respectively. In long- wave approximation the governing equations 
are found by a perturbation series in powers of the small parameter e^. We write in two 
dimensions (2D) 

Ui = + eui^ + e^Ui^ + . . . (la) 
Wi = Wig + ewi^ + e^Wi^ + . . . (lb) 
Pi = Pi, + tPi, + e^Pi, + . . . (Ic) 

where Ui and Wi stand for the x- and z-components of the velocities, respectively. The small 
parameter e = ^ 1 reflects the fact that the interface deflections are long scale, i.e. the 
mean film thickness is small compared to the typical lateral length scale A. 
Taking the large difference in length scales into account, it is natural to scale the system 
lengths like z = Hqz' and x = ^x', the velocities like Ui = uqu[ and Wi = euow'^, the time 
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like t = ^f, the pressures like Pi = the body forces like $i = and the 

normal and tangential interface forces like 11 = ^^^H' and r = ^j^t'. The primes denote 
the dimensionless variables, uq is a reference velocity of fluid 1 parallel to the substrate. 
Starting from the two-dimensional incompressible Navier-Stokes equations and the continu- 
ity equations for both liquid layers we derive the dimensionless equations in zeroth order in 
e. After substituting Eqs. in the governing equations and neglecting all terms of 0(e) or 
smaller, we drop the primes and the subscript zero and obtain for the scaled quantities in 
the lower layer, < z < h{x, t) 

dlu^ = d,P, (2a) 
dA = (2b) 
dxUi + dzWi = (2c) 

and in the upper layer, h{x, t) < z < d 

^ldlu2 = d,P2 (2d) 

d.P2 = (2e) 
dxU2 + dzW2 = 0. (2f) 



The variables 



Pi = Pi + <l>i and P2 = P2 + $2 (3) 



stand for reduced pressures which are the sum of the hydrostatic pressure Pj and the potential 
of the conservative body force $j (e.g. gravity force). The parameter /i = ^2/ ^^i represents 
the ratio of the dynamic viscosities. The boundary conditions at the lower and the upper 
boundaries read 

Ml = 0, Wi = at z = (4a) 



and 

U2 = 0, W2 = at z = d. (4b) 

The resulting interface conditions in zeroth order at z = h{x, t) are 

P-^ - P2 = AT + $ (5a) 

dzUi - fidzU2 = T (5b) 

dfh + Uidxh = wi (5c) 

dth + U2dxh = W2 (5d) 

ui = U2 (5e) 

Wi = W2 (5f) 

In the remainder of the paper we use the abbreviations A/" for the normal forces Eq. ()5ap. 
T for the tangential forces at the interface Eq. ()5b|) and $ = $1 — $2 for the body force 
potential. We mention that the occurrence of the body force potential $ in Eq. (j5aj) is caused 
by using reduced pressures Eq. (jSj). 

Because Pi and P2 do not depend on z (Eqs. ()2b|) and (j^ ) one can integrate Eqs. (Pa|) and 
(j2H]) in z. With the boundary conditions Eq. (jH) and the interface conditions Eqs. ()5b|) . 
the x-components of the velocities take the explicit form 

ui{x, z, t) = ^d^Pi + {-hd^M - hd^^ + B + T)z (6a) 
U2{x, z, t) = ^{d.Pi - dM - 9.$) [z^ - d^) + {z - d) (6b) 

o 1 -{d,Pi - 2dM - 29,$)/i/i2 + (a,Pi - DM - 5.$)(/i2 - rf2) - 2/ir/i 

Wltn iD = ; — ; 

2 {^-l)h + d 

where we used Eq. fjHajl to express dxP2 as a function of dxPi- 

Next, we derive the evolution equation for the interface profile h{x, t) and an explicit formula 
for the pressure d^Pi- To do so the continuity equation Eq. is integrated in 2;-direction. 
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With the interface condition Eq. (|^ and the chain rule we find 



ph{x,t) 

dth + dx / z, t) = 0. (7) 

Jo 

A similar evolution equation for h{x, t) is also derived for the upper layer using Eqs. pTj) and 
(l5d]) . Then Eq. (jS^) yields the identity 

/ rHx,t) r-d \ 

dx\ I ui{x,z,t)dz+ / U2{x^z,t)dz\ =0. (8) 

Jh{x,t) J 

To obtain Pi, Eq. (jHJ is integrated in x setting the integration constant to zero without loss 
of generality. This can be done if there is no additional lateral driving force as, for instance, 
in an inclined system. The resulting equation is solved using Eq. (^. The resulting pressure 
gradient is 

d^Pi = F,{h)d, (AT + $) + F2{h)T (9) 

with 

F^{h) = ^{d-hf{hii{Ad-h) + {d-hf) (10a) 
m) = ^id-h) (10b) 
D = {d-hY + hfi{h%fi - 2) + Adh^ - 6d^h + 4rf3). (lOc) 

The evolution equation Eq. ((Tj) can be written (with Eq. (|Uaj)) as 

dth = [Qi{h) (AT + $) + Q2{h) T] . (11) 

Using the same procedure we derive a similar equation for three dimensions 



a/i = V ■ 



Q,{h)V{M+^)+Q2{h)f (12) 



where V = {d^, dy). We note that in 3D an additional field occurs which can not be expressed 
as a function of h analytically. The importance and the equation of this mean fiow field are 
discussed in the appendix. However, we neglect the mean fiow in the following. 
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The mobilities are 

Qi{h) = \^ ' [d + /^(/x - 1)] (13a) 
Q,(h) = ^'^^ [h'i^i - 1) - d{d - 2h)]. (13b) 

The mobihty Qi{h) is positive for all values fi > and d > 0. It vanishes for h = and 
h = d. However, Q2{h) always changes its sign at 

he = — (14) 

The mobility Q2{h) does only exist in a system with shear-stress. In Fig.|21the mobilities 
are plotted for the parameters of Tab. H] with d = 1.3. There the change of sign of Q2{h) 
occurs at he ~ 0.91. 

In an one-layer system the viscosity of the upper gas layer is neglected. Therefore, taking 
the limit — > of the mobilities Eqs. ()13|) give the correct one-layer mobilities 

Qi(/^)iim = and Q2(/i)iim = (15) 

and the evolution equation Eq. ()12|) takes the well known form of the thin film equation for 
a single layer— 



dth = -V 
The limit of the pressure is 



— -h^^Pl lim + -h'^%im 



(16) 



A lim = Pi + $1 (17) 

since V-P2 = (see Eq. (|2H|) ) and $2 can be neglected as can be seen, for example, for gravity 
forces where $j oc pi and p2 <^ pi leads directly to $1 ^ $2- The same limit is reached by 
increasing the system thickness d — > 00. 
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B. Specific effects 



1. Gravitation and surface tension 

Incorporation of gravitation and surface tension (capillarity) provides the body force poten- 
tial and the normal interface force 

= (1 - p)Gh (18a) 
and M = -C-^V^/i, (18b) 

respectively. Thereby, G — {pighD/diiUo) is the gravity number, p — P2/ pi denotes the 
ratio of densities, — e^a/^iiiUo) denotes the capillary number and a is the dimensional 
liquid-liquid surface tension. 

2. Thermocapillarity 

The equations for the non-dimensional temperatures ©j (energy equations) read in zeroth 
order in e 

dlQi = (19a) 
9^02 = 0. (19b) 

The non-dimensional temperature is defined by 

T- -T 

e, = (20) 

where and Ti refer to the temperature at the upper and lower plate, respectively. The 
boundary conditions at the lower and upper rigid plate are 

Gi = 1 at z = (21a) 

and 02 = at z ^ d (21b) 
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and the interface conditions at z = h{x, t) read 

61 = 62 (22a) 
9,61 = A9,e2, (22b) 

where A = A2/A1 is the ratio of the thermal conductivities. Eqs. (fTII|l together with Eqs. (^1)1 
and give the temperature fields 

Assuming an arbitrary dependence of the surface tension a on temperature, the tangential 
interface condition Eq. ()5b|) has the form T = VS where S = ecr/ (/UiMq) is the dimensionless 
surface tension. Evaluation of VS at the position 2; = h{x,y,t) gives VS = dQT,-dhQ{h)-'Vh. 
If the surface tension depends linearly on temperature one gets 



r = M— ^^ Vh (24) 



where M = {—dTC AT e) / {^iUq) is the Marangoni number, 9^0" is the change of surface 
tension with temperature and AT = Ti ~ is the applied temperature difference. If the 
system is heated from below, M is positive for most fluids (normal thermocapillary effect). 
To derive the one-layer limit of the thermocapillary force Tiim one replaces A with a (con- 
ductive) Biot number A = Bi{d — 1) in Eq. (j24j) and takes the limit d —>■ 00. Considering 
the gas layer as a semi-infinite layer, we get the usual one-layer expression^ 

Bi 

^im = M Vh, (25) 

[Bi n + 1) 

where Bi 1. We note that the 'two-layer Biot number', Bi2, introduced by VanHook 
ET AL.— can be obtained by replacing A = ^^f^^'l~^^ in Eq. ()24|) and taking the limit d ^ 1. 
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This limit is valid since all dependencies on d are considered to be in Bi2 and therefore the 
neglect of the upper layer leads to d ^ 1. 



3. Disjoining pressure 

To avoid film rupture at the two bounding plates we incorporate disjoining pressures in the 
body force potential to model repelling stabilizing van der Waals forcesiS*^ 

Hi 



= ^ 



Z' 



^2 

= 



(26a) 

h 

(26b) 

z=h 



{d-zf 

The parameters Hi and H2 are Hamaker constants representing the interaction of the lower 
plate with liquid 2 through liquid 1 and of the upper plate with liquid 1 through liquid 
2, respectively^. Thereby, we neglect a part of the forces between the lower (upper) fluid 
and the upper (lower) substrate resulting from the finite thickness of the respective layer. 
The Hamaker constants determine the macroscopic contact angles. Since large macroscopic 
contact angles violate the used long-wave approximation we use Hamaker constants which 
provide a small contact angle. These corresponds to fluids which partially wet the substrates. 
Note, that for a very small distance between the plates, i.e. if both layers are ultrathin with 
thicknesses below 100 nm the disjoining pressure has the most important influence and all 
forces have to be included. Such systems are not the scope of the present work, but see^^ 
for a related system. 

4- Electrostatic field 

An electric field applied in z-direction is another way to cause structure formation. Consider 
two dielectric fluids with permittivities Si and £2, respectively. The upper and lower plates 
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serve as electrodes and a voltage U is applied. The vertical components of the electric fields 
in fiuid 1 and 2 read then in zeroth order lubrication approximation 

e2U 



El 

and E2 



e2h + ei{d — h) 
eiU 



(27a) 
(27b) 



e2h + ei{d — h) 

Horizontal components can be neglected at this order. Using the electrohydrodynamic stress 
tensor22: for pi =const., provides the effective electrostatic pressure at the interface by pro- 
jecting the stress tensor two times on the normal vector 



Pel = 2^o(^2 - ei)EiE2. 



(28) 



Scaling the voltage like U = V ^ piu^hQ/eoEie and dropping the primes, leads in zeroth 
order lubrication approximation to 

1 e{e - l)[/2 



el 



(29) 



2{eh + {d-h)Y' 

where Sq denotes the permittivity of vacuum and e = £2/siis the ratio of permittivities. In 
zeroth order in e shear stresses Tei are not present. 



C. Energy 



As for one-layer systems^ also here it is possible to express the r.h.s. of the evolution equation 
Eq. ()12|1 in variational form 



a/i = V • 



(30) 



corresponding to the evolution equation of a conserved order parameter field in a relaxational 
situation. Incorporating the above mentioned effects the energy E that corresponds to a 
Lyapunov functional can be written as 

"1 



E 



dx dy 



-C-\Vhf + V{h) 



(31) 
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with 



V{h) = ^(1 - p)Gh' + ^h'' + ^{d - hr' + Eth{h) + E,i{h), (32) 



and 

3M 



Ethih) 



A2(/x - Xfh In (/i) + (/i - A)2(d - h) In (d - h) 



2c/A(/i- A)2 
+A2(/i - -l) + d) In - 1) + rf) 



+ (A^/i - 2AV/i + AV(2/i -d) + 2Xfi{d - h) - f/{d - h)) In (/i(A - 1) + rf) 

(33) 

for thermocapillarity and 
for electrostatic fields. 

It can easily be shown2*24 that the Lyapunov functional E is monotonously decreasing in 
time {-^E < 0), if the mobility Qi{h) > which is always fulfilled. Note, that the free 
energy density for thermocapillarity Eth is a function of the ratios of viscosities /i, thermal 
conductivities A and layer thicknesses d. The dependence of Et^ on the viscosities, i.e. on a 
dynamical aspect of the system, does not occur in a one-layer system. 

III. RESULTS AND DISCUSSION 
A. Material parameters 

For our numerical investigations we focus on one specific two-layer system to allow for a direct 
comparison to experiments. We chose an oil-oil system used in Ref.-^, namely silicon oil 5cS 
and HT70. The parameter values are given in Tab.lH Note, that the given permittivities are 
only a rough estimate. 
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B. Linear stability 

To solve the hnear problem the normal mode ansatz 

h{x, y, t) = hk exp {ik^x + iky-y + xt)- (35) 
is used in Eq. (jl2p . Linearization provides the growth rate x 

X = -9Me{e-ki) (36) 



with ki = C 



e{e-iyu^ 3H2 Q2(l) ArfM 



(p - l)G + — - 3H, 



(37) 



is-l + df (d-l)4 Q,(l)^x-l + d)\ 

where kc is the cut-off or critical wavenumber and k"^ = k"^ + ky. The system is unstable 
for positive growth rates x > 0, i.e. for k < kc- Onset of the instability occurs with infinite 
wave length when kc = 0. 

1. Gravitation, surface tension and thermocapillary effects 

First, we study the situation without electric field, i.e. U = 0. For Hi,H2 ^ 1 the linear 
stability is determined by p and M only. In the isothermal case (M = 0) the system is 
unstable for p > 1, i.e. the system is gravitationally or Rayleigh- Taylor (RT) unstable. In 
the heated case, Eq. ()36|) provides the critical Marangoni number 



2(A- l + rf)2(d- l)(p + d- 1) / 3H._ 



We note that the critical Marangoni number correspond to the one found by Smitb^^ for thin 
layers of viscous liquids. Inspection of Eq. (|36|) shows that the sign of the thermocapillarity 
term does not only depend on the sign of the Marangoni number but also on the sign of 
the mobility ^2(1). This implies that the sign of ^2(1) determines whether M must be 
larger or smaller than M^. to get an instability. Denoting the zero crossing of Q2{h) by he 
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(see Eq. fll4|) ). one finds that for he > 1 the Marangoni number M lias to be increased over 
Mc for tlie system to become unstable, whereas for he < 1 it has to be decreased below 
Mc- The destabilizing direction of heating is determined by substituting /i^ = 1 in Eq. (fT^. 
Instability results if 

fx < {d- if for M> Me (39a) 

and if fi > {d - if for M < M^. (39b) 

The critical system thickness (critical viscosity ratio, respectively) reads then 

4 = v^+l ifie=id-lf) (40) 

as already found in Ref.— . The dependence of the growth rate on the wavenumber is shown 
in Fig.Elfor three different situations at a system thickness d = 1.3, parameters from the first 
column in Tab.U and Hi = H2 = 0. Because p = P2/P1 > 1 without heating the system is 
Rayleigh- Taylor unstable (solid line). Since the system thickness is smaller than the critical 
one de = 1.43 (Eq. ()4()|l ). heating from below with M > Me = 0.89 damps the RT instability 
(dashed line). As indicated in Eq. ()39b|) heating from above amplifies the Rayleigh- Taylor 
instability (dotted line). The stabilization mechanism for M > Me is directly correlated with 
the sign of the mobility Q2{h). The zero crossing of Q2{h) is at he ~ 0.91 for the chosen 
material parameters. Therefore, the mobility Q2 is positive in the linear regime {h ^ 1) and 
heating from below (M > 0) acts stabilizing. 

The stability diagrams in Fig. ^ show the critical Marangoni number Me (Eq. ()38|)) in de- 
pendence of the ratio of viscosities fi for d = 1.3 and A from the first column in Tab.Q] 
{Hi = H2 = 0) . The left (right) panel represents a system that is Rayleigh- Taylor unstable 
(stable) for M = 0. At the critical viscosity fie = 0.09 the critical Marangoni number Mc 
changes its algebraic sign in both cases. Obviously, thermocapillarity dominates for strong 
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heating (large \M\), i.e. RT is negligible. 

To understand the mechanism of the stabilization of a RT instability for M > and /i > /Xc, 
we consider a small deformation of the interface in negative 2;-direction as sketched in Fig.El 
First, we discuss the isothermal case (M = 0) where for p > 1 the system evolves due to 
its Rayleigh- Taylor instability. The viscous time scales of the two layers are responsible for 
the direction of the fluid velocity in the vicinity of the deformation minimum. For d > dc 
the viscous time scale of the lower layer ti oc h^/fii is faster than the viscous time scale 
T2 oc {d — h^Y / ^2 of the upper one. Therefore, the lower layer is the driving layer and 
velocities are directed away from the deformation minimum (solid arrows in Fig.|Sl(a)). For 
d < dc the velocities are directed towards the deformation (solid arrows in Fig.|Sl(b)). 
When heated from below (M > 0) the temperature is highest at the deformation minimum. 
The accompanying surface tension gradient causes a flow away from the minimum (dashed 
arrows in Fig.|SJ. This leads to an amplification of the perturbation for d > dc, because 
thermocapillarity acts in the same direction as the Rayleigh- Taylor mechanism (Fig.El(a)). 
For d < dc thermocapillary forces act in the same direction as before, but the driving of 
the upper layer leads to flow towards the deformation minimum (Fig.|31(b)). Therefore 
thermocapillarity damps out the Rayleigh- Taylor instability. 

Finally, we display in Fig.|Hl the critical Marangoni number Mc in its dependence on the 
system thickness d as calculated from Eq. ()38|) . To avoid a destabilizing influence of a 
Rayleigh- Taylor instability we take the parameters from the second column in Tab. HI i.e. we 
interchange the two liquids. For d > dc ^ 3.34 a minimum in Mc is observed. It reflects 
the antagonistic influences of the system thickness and the temperature gradient in the 
lower layer {driving layer). For large d the temperature gradient in the lower layer tends to 
zero (from Eqs. ()21|) and (f^Hj) ). implying a large critical Marangoni number. Decreasing the 
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system thickness leads to decreasing M^. When decreasing d further the mobihty Q2{h = 1) 
tends to zero. Hence, for d ^ one again finds an increasing Mc. 

For system thicknesses d < d^ one observes a monotonously decreasing |Mc| for d ^ 1. In- 
cluding disjoining pressures as repelling forces changes the behavior for very small thickness 
of the upper layer d — 1 qualitatively, but has no influence otherwise (dashed line). Specifi- 
cally, for d < dc the Hamaker constant of the upper layer H2 (representing the interaction of 
the upper substrate with liquid 1 through liquid 2) forces an extremum of Mc. Decreasing 
the thickness of the upper layer towards d ^ 1 the stabilizing Van der Waals interaction 
finally dominates allowing to consider systems with very small d—1 as stable. The Hamaker 
constant Hi (representing the interaction of the lower substrate with liquid 2 through liquid 
1) causes only a slight shift of Mc. It does not change the extremum in this region. 

2. Electrostatic field 

We conclude the discussion of the linear stability by regarding the influence of a vertical 
electrical field only. Neglecting thermocapillarity (M = 0) and gravitation {G = 0) yields 
for Hi, H2 1 the critical voltage from Eq. ()36|) 



Note, that the direction of the applied voltage (±2;) has no influence on the stability. Using 
parameters of the second column in Tab. |l] with Hi = H2 = 0.01 and d = 4 we find the 
critical voltage Uc ~ 4.5. A voltage of [/ = 30 provides the wavenumber kmax ~ 0.18 for the 
maximal growth rate x- We use these parameters to study the time evolution with the fully 
nonlinear equation below in Section Fill D 31 




(41) 
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C. Implications of the variational formulation 

The variational formulation of the evolution equation (Eq. ()30|) ) provides an energy or Lya- 
punov functional based on a gradient energy and a free energy density V{h) (Eq. (j32j) ) . The 
latter can be used in a Maxwell construction that allows for the prediction of the character 
of the resulting structure in the long-time evolutioni^ as well as the study of metastable 
statesii. 

1. Maxwell construction 

First, we want to determine whether holes or drops are formed in the long-time evolution. 
Since a Lyapunov functional Eq. (|31|) exists, the final equilibrium thickness profile corre- 
sponds to its global minimum. The mean height ho is a conserved quantity, i.e. an increase 
of the interface height in any region is accompanied by a decrease somewhere else. When 
minimizing the energy functional, this constraint has to be taken into account by a Lagrange 
multiplier, A^,, namely by supplementing the free energy density V{h) by the term A^/;,. 
Assuming a very large system in a late stage of coarsening the gradient term of the energy 
functional can be neglected and the local free energy suffices to derive the long-time behavior. 
First, consider a V{h) possessing one minimum and a monotonously increasing slope, i.e. a 
V{h) with positive second derivative everywhere. Then, deforming the interface increases 
the local part of the free energy functional (which is further increased due to the gradient 
term) since the energy loss is due to mass conservation accompanied by a larger energy gain. 
Therefore, for such a V{h) the (energetically minimized) final solution is a flat interface. 
Contrary to this, a V{h) with two minima may allow to minimize the local free energy (and 
may even overcome the energy gain due to the gradient term) by deforming the interface 
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since the free energy may decrease for both h > ho and h < ho. In this case the flat interface 
is hnearly unstable and the system will realize two film thicknesses {hi < ho = 1 < h2). 
In analogy to spinodal decomposition the two film thicknesses can be seen as two different 
phases, and the evolution of the film thickness profile corresponds to a phase separation^. 
Mathematically formulated, the phase separation occurs if it is possible to find a double 
tangent, where the curve V{h) lies everywhere above this tangent. The slope of the tangent 
corresponds to the Lagrange multiplier A^, and the points where it touches the curve V{h) 
give the two equilibrium values of h. In the present case the existence of the double tangent 
is assured by the stabilizing disjoining pressures. Without the latter the equilibrium film 
thicknesses may be found outside the gap between the two substrate indicating finite 'true' 
contact angles. The double tangent condition is equivalent to a Maxwell construction in the 
[—dhV {h)]-[h] space. 

Resulting from mass conservation, the ratio of the surface areas S = S1/S2 of the two 
equilibrium film thicknesses {hi < ho and h2 > ho) defines the solution morphology. In 
accordance with observed structures in thin films we call S < 1, S > 1 and 5* = 1 hole, 
drop and maze solutions, respectively. However, for systems with 5* close to 1 but 5* < 1 
(5 > 1) the solution reveals its visible hole (drop) character not until the final stationary 
state. Note, that we use the expression 'hole' ('drop') for a hole (drop) in (of) the lower 
fluid. Obviously a hole (drop) in (of) the lower fluid corresponds to a drop (hole) of (in) the 
upper fluid. 

We focus on the system thickness d as control parameter for the phase separation since d 
can be controlled easily in experiments. Fig.[7|(a) shows the plot of —dhV{h) for different 
system thicknesses d. The Maxwell point hM (via a Maxwell construction) provides then 
a criterion for holes or drops. If hM > 1 (5* > 1) drops are preferred (dashed line), in the 
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other case (5* < 1) holes are expected (sohd hne). For d = 2 the transition from holes to 
drops takes place (dotted hne). We mention, that this Maxwell point does not coincide 
with the critical system thickness dc- 

2. Metastability 

As stated in the linear investigation above in Section llllB the Rayleigh- Taylor (RT) insta- 
bility can be stabilized by heating from below (for d < dc). This is also confirmed by a fully 
nonlinear integration in time (not shown). Nevertheless, a RT unstable system stabilized 
by thermocapillarity can be metastable, i.e. it may be nonlinearly unstable to (large) finite 
perturbations. This metastability can also be studied using a Maxwell construction as shown 
in Fig.[7|(b). Under isothermal conditions (dashed line) the system is RT unstable indicated 
in Fig.[7|(b) by the fact that the vertical line at /i = 1 crosses the dashed curve in-between 
the two extrema. When heating from below with M = 10 (solid line) the system is linearly 
stable. However, the Maxwell plot has still two extrema. Since the mean system thickness 
{Hq = 1) is situated to the right of the maximum and the local free energy (y{h = 1)) is 
larger then the free energy of the Maxwell point {V{hm)) the system is metastable. 
Fig. IHl shows the critical Marangoni number Mc in dependence of the system thickness d. For 
d < dc {d > dc) heating from below (above) with M > Mc (M < Mc) stabilizes the system. 
However, a Maxwell construction shows a metastable state for all the displayed Marangoni 
numbers (|M| < 60). This metastability was also found in experiments with oil-air layers^i. 
A qualitative understanding of the metastability is given by the zero crossing he of the 
mobility Q2 (Eq. ()14|) ). For perturbations larger than |1 — hd the interface is destabilized 
since both gravity and thermocapillarity destabilize the system. 

Fig. El shows a snapshot from a two-dimensional numerical run for a RT instability without 
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thermocapillarity (dashed line). Heating from below (M = 10) stabilizes the system and 
leads to a flat stable interface for small perturbations (not shown). However, starting with a 
finite perturbation in the vicinity of x = 50 leads to a state with a local strong modulation 
(solid line). 

D. Long time evolution 

Three-dimensional numerical integrations of the nonlinear equation Eq. ()12|) are done with 
an ADI method (Alternating Discretization Integration). In the first half time step the lin- 
ear part is integrated implicitly in x-direction, in the second half time step in ^/-direction. 
The nonlinear part is calculated explicitly. We use periodic boundary conditions in x and 
y and initially disturb the fiat interface with small random fluctuations T]{x,y,t). Thereby 
the average height is conserved (/ 1 + r]{x,y,t) dxdy = Hq = 1). Further on, we distinguish 
between short-time evolution and long-time evolution. Roughly speaking linear effects de- 
termine the dominant length scales of the short-time evolution. Nonlinear effects dominate 
the long-time evolution that is characterized by coarsening processes. 

1. Ray leigh- Taylor instability 

In the isothermal case without electric fields, gravity is the only possibly destabilizing influ- 
ence. The long-time evolution of a system with d = 3 > dc, G = 5 and material parameters 
from the first column in Tab.lJ is shown in Fig.Qni Initially small disturbances of the fiat 
interface evolve into a drop structure {t ^ 1000). For larger times small drops vanish due to 
coarsening and finally the system settles at the global energetic minimum corresponding to 
one large drop (not shown). Fig.lTTl gives the evolution for d = 1.3 < dc and G = 20. Here, 
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the short-time evolution results in a hole pattern (t < 200). Subsequently, the long-time 
coarsening towards structures of larger extension sets in {t = 1000) and finally ends with 
one large hole (t > 2 x 10^, not shown). 

The use of G = 20 for d = 1.3 and not G = 5 as for (i = 3 assures a smooth growth in the 
short-time evolution. By 'smooth' we mean a gradual growth of all holes until they have 
rather large amplitudes. Using instead G = 5 for = 1.3 results in a rapid non-smooth hole 
evolution in-between the short-time and long-time domain, i.e. the structure is determined 
by the linear wavelength at the very beginning of the evolution only. As soon as nonlinear 
terms become important {\ri{x,y.t)\ <^ 1 is violated) only part of the linearly developed 
structure evolves. Here, this rapid hole evolution is caused by the stabilizing mechanism 
of the disjoining pressure at the upper plate (Hamaker constant H2 = 0.01) which can no 
longer be neglected even for small perturbations of the fiat interface. This affects both the 
linear and the nonlinear evolution of the interface. Namely, it suppresses interface evolutions 
for h > 1 and therefore causes a rapid evolution for /i < 1. We note again, that here the 
rapid hole evolution is caused by disjoining pressures, whereas the mobility Qi{h) has no 
effect. 

Fig. ^1 displays the time evolution for d = 2, G = 5, and parameters from the first column 
of Tab.lH For these parameters, neither drops nor holes are energetically preferred. This 
leads to a clearly visible maze or labyrinth structure which also shows the typical coarsening 
dynamics at long times. All long-time solutions shown (drops for d = 3, holes for d = 1.3 
and maze structures for d = 2) correspond to the predictions of the Maxwell construction 
in Fig.[7| 
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To quantify the coarsening behavior we calculate at each timestep the mean wave number 



tors are distributed on a small annulus. Therefore the approximation (A;)^ ^ (fc^) holds and 
the mean wave number (k) can also be taken as a qualitative measure of the mean curvature 
of h{x,y,t). Fig.lT^ shows the dependence of {k) and of the corresponding energy (Eq. (jSH)) 
on time for the three numerical evolutions discussed above. The energy decreases always 
monotonously in time as expected. The mean wave number (k) shows two local extrema at 
tjnin and tmax, respectively, with t^m < tmax- In the region between the two extrema the 
amplitude of the interface deflections outgrows the linear regime. Further on, the averaging 
in Eq. ()42|) allows to interpret the strength of the maximum as a measurement for global 
amplitude growth. The absence of a local maximum indicates that all linearly evolved drops 
or holes evolve globally and uniformly towards larger amplitudes, whereas a strong peak 
indicates a nonlinear evolution of a few linearly evolved drops (holes) only. 
The mean wave number at the local minimum corresponds to the wave number of the max- 
imal growth rate from the linear investigation (thin horizontal lines). Hence the evolution 
for if: < train IS determined by linear terms, i.e. the wavenumber with the maximal linear 
growth rate emerges in the system. Therefore, the region around the two extrema can be 
considered as the frontier between short-time and long-time evolution. In the long-time 
evolution (t > tmax) nonlinear effects dominate and a scaling lawi^ 



can be extracted which reflects the coarsening of the system for long times. To determine 
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(42) 



where h{kx, ky, t) denotes the Fourier transform of h{x, y, t). We mention, that the wavevec- 



{k) = c- 



(43) 



a 'true' scaling exponent a statistically significant average of many numercial runs with 
different initial conditions is necessary. Due to the strongly time-consuming character of 
the necessary computer calculations only a few runs were used to determine the respective 
tendencies of scaling exponents presented here. However, in the following we call them 
shortly 'scaling exponents'. 

We find a nearly identical scaling exponent of /5 ~ 0.14 for drops {d = 3, solid line, tmax ~ 
150), holes {d = 1.3, dashed line, tmax ~ 210) and maze structures {d = 2, dotted line, 
tmax ~ 320). Neither the system thickness d nor the gravity number affect the long-time 
scaling. 

To identify the acting coarsening modes we illustrate the flow pattern by calculating differ- 
ences between the interfaces h{x, y, ti) and h{x, y, at different times ti and ^2, respectively, 
shown for the evolution of drops in Fig.El (corresponding to Fig.llO|). One can identify two 
different coarsening mechanisms being dominant at different times within a single long-time 
evolution. As a result of the short-time evolution many small drops exist. Neighboring 
drops attract each other strong enough to move the entire (small) drops and finally combine 
to one large drop sitting at an intermediate position. This translational coarsening mode 
is illustrated in Fig.El(a) where its signature in the difference plot is that all drops have 
white (mass gain) and black (mass loss) parts of their edges. 

For larger times a transition from the dominant translational mode to a dominance of the 
volume transfer mode takes places. Now the mean distance of the drops is too large to get 
the (large) drops moving. Only material is transported between the sitting drops resulting 
in a slow disappearance of smaller drops and the growth of the larger drops. This mass 
transfer mode is illustrated in Fig.El(b) where its signature in the difference plot is that 
there exist drops that have completely white or completely black edges. 
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Finally, we want to show that our descriptive explanation for the viscous timescales (and 
therefore for the interface velocity directions) given above in Section Fill Bl is in concordance 
with the fully nonlinear evolution. We integrated Eq. numerically in two dimensions 
for a RT unstable system with Hi = H2 = (integration was stopped before film rupture 
occurred). The x-components of the velocities Ui{x,z,t) and U2{x,z,t) can be calculated 
from Eq. ©. They are plotted in Fig. (right plots) for the position x = 20. The left panels 
of Fig.^1 display the interfaces and the contour lines of the streamfunction ip {u = —dzf, 
solid lines ip > 0, dashed lines (p < 0). For fi < fic (Fig.UHKa)) the interface velocity is 
negative, i.e. fluid moves away from the deformation minimum and therefore the lower layer 
is the driving layer (solid arrows in Fig.|Sl(a)). For yU > /ic (Fig.^l(b)) the interface velocity 
is positive. Hence fluid moves to the deformation minimum and the upper layer is the driving 
layer (see also solid arrows in Fig.|31(b)). 

2. Thermocapillary effects 

In this section we include thermocapillary effects (using Eq. (j2H)). Thermocapillarity can act 
both stabilizing and destabilizing as seen from linear analysis. To avoid an amplification due 
to gravity we use parameters from the second column in Tab.|ll Therefore, p < 1 and gravity 
stabilizes the fiat film. Eq. ()39b|l provides then the critical system thickness d^ ^ 3.34. An 
unstable initial fiat film is obtained for d > d^ {d < d^) by heating from below (above). In 
the following, we use d = 2 and d = 4 to illustrate the destabilization by different directions 
of heating. 

For d = 2, the critical Marangoni number is Mc ~ —5 and the system becomes unstable for 
M < Mc- Fig.lTHl shows a numerical run for M = —10. Small initial disturbances evolve 
smoothly into drops which coarse in the long-time evolution [t = 3 x 10^) corresponding to 
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the prediction of a Maxwell construction (Maxwell point Km ~ 1.34). Fig.^l (solid line) 
shows the mean wave number {k){t) and the energy versus time. The energy is again a 
monotonously decreasing function in time. The mean wave number {k){t) has again two 
local extrema. The minimum reached at tmm ~ 1800 corresponds to the fastest linear wave 
number. The long-time coarsening sets in after the local maximum at t^ax ~ 4000 and the 
scaling coefficient defined in Eq. ()43|) is determined to be /5 ~ 0.16. 

For d = 4 the critical Marangoni number is Mc ~ 37.84 and we use M = 70 for the numerical 
run displayed in Fig.^J Starting from small perturbations one hole evolves rapidly at 
t ~ 700. Subsequently, more and more holes arise {t = 1100). For t > 1100 coarsening sets 
in and in the long-time limit a single hole remains (t = 2 x 10^). This corresponds to the 
prediction of the Maxwell construction (Maxwell point /ia/ ~ 0.51). The mean wave number 
(k) in Fig.EI (dashed line) shows again a minimum corresponding to the fastest linear wave 
number and a very pronounced maximum indicating the rapid evolution of one hole between 
the two extrema. Since the averaging in Eq. ()42|) gives approximately also the root of the 
mean curvature the abrupt rise of (k) is obvious even for the evolution of a single hole. The 
long-time scaling is with (3 ~ 0.27 remarkably faster than for d = 2. Note, that we found 
numerically that the scaling exponent for d = 4 does nearly not depend on the Marangoni 
number. 

The differences in the short-time evolution for d = 2 and d = 4 (rapid evolution of one hole 
for d = 4 versus smooth evolution of many drops for d = 2) can be understood in terms of 
the effective mobilities Q[{h) = G{l-p) Qi{h) and Q'^{h) = MQ2{h). shown in Fig.m For 
d = 4, Q2 crosses zero close to h = 1 (thick lines, he ~ 1.2 from Eq. p4|) ) and Q'l increases 
for increasing h. Therefore, the interface evolution is slowed down for h > 1 (and finally 
stopped for h > he) and accelerated for h < 1. This results in a rapid hole evolution. For 
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d = 2, both mobilities show an approximate symmetry around h = 1 (thin hnes). Hence, no 
interface thickness is suppressed allowing for a smooth evolution of drops. 

3. Application of an electric field 

Finally, we illustrate the time evolution caused by a vertically applied electrical field. We 
use the parameters from the second column in Tab.|l] for an isothermal (M = 0) system 
with G = 0. Fig. 1201 shows snapshots of the long-time evolution for an applied voltage 
U = 30. Initially small disturbances of the interface evolve smoothly to drops and the 
long-time coarsening sets in at t ^ 10000. The dependence of the mean wave number on 
time shown as dotted line in Fig.^l shows a minimum at tmin ~ 1500 and a maximum at 
tmax ~ 4700. Again, the minimum coincidences to the wave number of the fastest linear 
mode kmax ^0.18. 

The derived long-time scaling exponent f3 ~ 0.04 is small compared to the exponents mea- 
sured above for the Rayleigh- Taylor and thermocapillary instabilities. In absolute values we 
find only a small change from (k) = 0.18 at t = 10^ to (k) = 0.16 at t = 10^. We conclude, 
that one can consider the length scale of the pattern in the long-time evolution as frozen to 
the value of the wavelength of the fastest linear mode {2tt/ kmax) at least up to t = 10^. 

IV. CONCLUSION 

Using long-wave approximation we have derived a single evolution equation for the interface 
profile of a two-layer system bounded by rigid plates. This equation is written in a general 
form to facilitate the inclusion of arbitrary body forces and normal or tangential forces at 
the interface. In the analysis of the model presented here, we have focused on the influences 
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of gravity, thermocapillarity and electrostatic fields. 

We have shown that the mobility of the normal-stress and the body force terms Qi{h) is 
always positive as in the one-layer case. However, it tends to zero not only for h ^ as 
for the one-layer system but also for h ^ d. The latter has no counterpart in one-layer 
systems where the mobility increases monotonously with the film thickness^. The second 
qualitative difference is the sign change of the mobility for the shear-stress term Q2{h). 
Both mobilities can affect strongly the linear and nonlinear evolution. For instance, the 
direction of heating needed for destabilization is determined by the zero crossing of the 
mobility Q2{h). Furthermore, the shapes, zeros and extrema of the mobilities allow at least 
a qualitative prediction of the dynamics of the system without any numerical investigation. 
Moreover, in contrast to weakly nonlinear theories^! we are able to check these descriptive 
criterions integrating the fully nonlinear equation numerically. A non-smooth rapid hole 
(drop) evolution in-between the short-time and long-time regime can already be estimated 
from the trend of the mobilities. 

Remarkably, although in the heated case the system is dissipating energy through convection 
within the drops (or around the holes) even when the final stable state is reached, the use of 
long-wave approximation allows for a variational formulation using an energy or Lyapunov 
functional for the film thickness profile. The film thickness evolution equation takes then 
the form of the simplest possible equation for the dynamics of a conserved order parameter 
j|g| j36^37^38 prominent representative of this class of systems is the Cahn-Hilliard equation 
describing the evolution of a concentration field for a binary mixture^. In contrast to the 
one-layer case^i*^, here the energy itself depends on material parameters that characterize 
the dynamics of the system, namely the ratio of viscosities. We used the energy functional 
to predict the expected long-time behavior, namely the evolution of holes, drops or maze 
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structures. It also allows for the study of metastable states. The predictions have been 
confirmed by fully nonlinear numerical integrations of the evolution equation. 
Using a linear stability analysis we have discussed the conditions for a gravitational or 
Rayleigh- Taylor instability, thermocapillary destabilization and stabilization, and an elec- 
trohydrodynamic instability for dielectric liquids under the influence of an electrical field. 
We have shown that thermocapillarity may act stabilizing as well as destabilizing depend- 
ing on material parameters. The behavior becomes intuitively clear because when treating 
both layers in the same way no direction of heating should be preferred. This implies that 
depending on material parameters both ways to destabilize the system - heating from above 
and heating from below - have to be possible. 

We have given special emphasis to the study of the possibility to stabilize a Rayleigh- 
Taylor unstable two-layer system by heating from below. This seemingly counter-intuitive 
behavior first discussed in Ref.— is a typical property of multi-layer systems and is directly 
correlated with the change of sign of the mobility Q2{h). However, we have shown that 
the stabilized Rayleigh- Taylor system is metastable. This explains a problem encountered 
in the experiments of BuRGESS et al.^^. Although they could stabilize a Rayleigh- Taylor 
unstable 'oil on air' system by heating from below this was only possible in ten percent 
of the experimental runs. This is due to the fact that the preparation of the initial flat 
film involved large amplitude disturbances. Because of the metastability of the system this 
results in the destabilization of ninety percent of the runs because the disturbance is larger 
than the critical one. The transition from the two- to the three-dimensional equation has 
shown that a weak mean flow arises. We neglect this additional flow in the main part of the 
present work due to its very weak influence on the presented results. However, the mean 
flow has been discussed in detail in the appendix. 
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We have implemented numerical schemes for two- and three-dimensional versions of the 
fully nonlinear equation and have given an overview showing different possible long-time 
evolutions consisting of coarsening hole, drop or maze patterns. Furthermore, we have 
analysed the scaling behavior by calculating the time dependence of the mean wave number 
of the patterns and extraction of the tendency of scaling exponents. 

Isothermal two-layer systems, i.e. taking into account gravitation and disjoining pressures 
only, show the same long-time scaling for different system thicknesses d. However, incor- 
porating thermocapillarity the system thickness d affects the long-time scaling essentially. 
This scaling behavior is in contrast to the one-layer case which was found to be determined 
by one scaling exponent {j3 ~ 0.21)^°. The one-layer coefficient lies within the two-layer 
range 0.16 < /5 < 0.27 found here. Moreover, an isothermal (M = 0) and agravic {G = 0) 
electrohydrodynamically unstable system reveals a very small scaling coefficient {f3 ~ 0.04). 
To our knowledge for this class of evolution equations such a slow long-time scaling is found 
for the first time. 

Finally, we have shown that tangential interface forces (thermocapillary forces) allow for 
rapid hole (or drop) formation after the short-time evolution. Again, this mechanism can 
be understood in terms of the change of sign of the tangential mobility Q2{h). 
In the present work, disjoining pressures were solely used to inhibit rupture of the layers 
to allow for the study of the long-time behavior. This corresponds to the assumption that 
liquid 1 and 2 wet the lower and upper plate, respectively. However, our present theory is 
not apt to describe situations where both layers are ultrathin (less than 100 nm), a situation 
gaining more and more importance for research communities and industrial applications. 
Then disjoining pressures dominate and the used terms are not exact enough because part 
of the forces between the lower (upper) fluid and the upper (lower) substrate are neglected. 
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Furthermore, the used Van der Waals interaction should be supplemented by additional 
short-range interactions^^. A further analysis of disjoining pressures in two-layer systems 
seems promising and will be published elsewhere. 

We emphasize our results for the action of a vertical electric field since recent experiments 
focus on such systems, as for example, done by LiN et al. using two polymeric liquidsSSiSi. 
They monitored the time evolution up to the impingement of the lower polymer layer on 
the upper electrode and showed a series of snapshots of the evolving morphology (Fig. 4 of 
Ref.— ). Interestingly, they found a nearly constant length scale of the evolving columnar 
structures from the early stage on, corresponding to the fastest linear mode. This corre- 
sponds to the results of our linear and nonlinear analysis of this case that we performed for a 
comparable ratio of permittivities. Furthermore, a visible concordance of both timeseries ex- 
ists. The small scaling exponent (3 ~ 0.04 found here can be regarded as a structure length 
frozen to the fastest linear wavelength. This indicates, that our model gives reasonable 
results even for macromolecular liquids. 

Finally, we stress the advantages for physics as well as industrial applications of bounded 
two-layer systems. From a physical point of view our single interface equation captures 
both the long-wave evolution and the interface interactions of two fluids. Therefore it allows 
for detailed investigations how the fluid properties of both the upper and the lower fluid 
layer determine the stability, metastability and short-time as well as long-time evolution. 
Furthermore ultrathin films play already a major role to create desired structures or stable 
flat interfaces. Usually one-layer equations are used for that kind of industrial applications. 
However, controlled boundary conditions, well defined bulk properties and consequently well 
defined interface actions enhance the accuracy of experiments as well as the examination of 
theory and experiment. 
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APPENDIX: THREE-DIMENSIONAL EVOLUTION EQUATION 

The generalization of the derivation of the film thickness equation to three dimensions is 
straightforward up to Eq. (jS} for which the three-dimensional (already integrated) version 
reads 



where V = {dx,dy) and 



A{h) (yPi - Fi(/i)V(Ar+ $) - F2(/i)f)] = (A.l) 



Aih) = ^ ^ (A.2) 



with D from Eq. (fTUc|) . 

Eq. ()A.1|) cannot be solved analytically for the pressure Pi. However, to fulfill Eq. (TOjl its 
argument must be of the form 

^ Yoi{fe,) = VPi-Fi{h)V{U+^)-F2{h)f (A.3) 



m 

with / as an arbitrary scalar function, rot(/e^) = {dyf, —dxf,0) and Fi{h) from Eq. pOap . 
Note that this function / is already present without symmetry breaking effects (e.g. inclined 
or rotating systems). Substitution of the gradient pressure from Eq. ()A.3j) in the three 
dimensional correspondent to Eq.© provides the evolution equation for h 



d,h = V ■ 



Qi{h)V{M +^)+Q2{h)f +V ■[Q?Xh)roi{fe,)] (A.4a) 



where / as a function of h and its spatial derivatives is given implicitly by 

--^A/ = dh(^-^^ rot(/ie,)-rot(/e,) + 9,Fi(/i)rot(/ie-;)-V(Ar+<|.) 

+dhF2{h) vot{h e,) ■ f + F^ih) rot(f) ■ e,, (A.4b) 

resulting from applying the curl to Eq. ()A.3|) . 
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The third mobihty 



Qsih) = Fi{h) - 1 



(A.5) 



is a monotonically decreasing function in h (Qsi^O) = 0,Q3{d) = —1). Note that dhQsi^h) = 
dhFiiJi) is zero for h = and h = d and negative for < h < d. The additional function 
/ reflects the possible mean flow of the system induced by a vertical vorticity contribution. 
We mention that in the one-layer limit (/i ^ or c/ — > oo) this vorticity contribution is zero 
{1/A{h) — > 0) and the usual normal condition for the pressure is recovered. 
Substituting Af, $ and T (from subsection I11B|) in the rhs of Eq. ()A.4b|) provides 



showing that solely the surface tension contributes to /. 

Due to the properties of dhFi{h) the vorticity / is basically determined by the geometrical 
shape of the interface. It can be shown with Eq. ()A.6|) that rotationally symmetric structures 
are not affected in the strongly nonlinear regime by this additional vorticity contribution. 
Furthermore Eq. ()A.6|) reveals that for small deviations from symmetric states the vorticity 
/ acts in the same direction as surface tension itself in Eq. ()A.4a|) . 

We have checked and confirmed all numerical runs of this paper taking into account Eq. ()A.6|) . 
In general the contribution of / is small. Therefore we neglect it and consider the evolution 
equation (fT^ only. However, we note that the effects of Eq. ()A.6|) on the evolution equation 
are an attractive and important subject for future investigations. 



A{h) 



1 




rot{he,) ■ rot(/e;) = C-^dhFi{h) rot{he,) ■ V(A/i) 



(A.6) 
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r mid Z - r luici i 


til lu - Silicon oil bcb 


Silicon oli oco - til (\J 


density p = P2/P1 


1 OOi^ 
i.OZD 


U.04o 


viscosity /i = 


U.ioZD 


/I Q 


thermal conductivity A = A2/A1 


0.598 


1.67 


permittivity e = £2/^1 


0.77 


1.3 



TABLE I: Material parameters for a silicon oil 5cS - HT70 system taken from Ref. . The values 
in the first column are for HT70 on silicon oil 5cS, whereas for the second column the liquids are 
interchanged, i.e. silicon oil 5cS on HT70. 
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Figure captions: 

Fig.^ Sketch of the system. The two layered immiscible liquids are bounded by two rigid 
smooth plates. The flat interface is situated at the mean height ho- The position 
h{x, y, t) of an evolving interface profile is a function of x, y and t only. 

Fig. 121 Shown are the mobilities Qi{h) (normal-stress and body force terms) and Q2{h) 
(tangential-stress term) for (i = 1.3 and /i = 0.1826. The mobility Qi{h) is always 
positive and vanishes for /i — >• and h ^ d. The mobility Q2ih) changes sign at 
he ^ 0.91. 

Fig. 121 The dependence of the linear growth rate x oii wavenumber k ioi d = 1.3, p = 1.826, 
fi = 0.1826 and A = 0.598 (parameters from the first column in Tab.jH) at G = 1 and 
= 20 {Hi = H2 = 0). In the non-isothermal case M = ±1. 

Fig. 01 The stability diagrams show the critical Marangoni number in dependence of the 
ratio of viscosities /i for = 1.3, A = 0.598 and G = 5 (Hi = H2 = 0). > 
(Mc < 0) corresponds to heating from below (above). The direction of heating leading 
to destabilization changes at the critical viscosity pc = 0.09. The left (right) panel 
corresponds for M = to a Rayleigh- Taylor unstable (stable) system with p = 1.826 > 
1 (p = 0.548 < 1). 

Fig. El Sketch illustrating the mechanism of the (a) destabilizing and (b) stabilizing thermo- 
capillary action for a Rayleigh- Taylor (RT) unstable system heated from below. Solid 
arrows indicate the liquid velocity close to the interface in the respective driving layer 
for pure RT, dashed arrows signal the effect of thermocapillarity. (a) Due to thermo- 
capillarity the total flow in the driving layer is directed away from the deformation as 
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known from one-layer systems. This causes an amplification of the disturbance, (b) 
The viscous timescale oc {d — /io)^/At2 of the upper layer is faster than the one of 
the lower layer. The influence of the upper driving layer dominates leading to flow to 
the deformation minimum. Thermocapillarity causes flow in the opposite direction, 
thereby weakening or completely damping the disturbance. 

Fig. El Shown is the critical Marangoni number Mc versus the system thickness dfoi p = 0.548, 
/i = 5.48, A = 1.67 and G = 5. The critical system thickness is dc ~ 3.34. For d < dc 
heating from above acts destabilizing and is a monotonic function of d. For d > d^ 
heating from below destabilizes and a minimum results from competing mechanisms 
(see main text). Inclusion of stabilizing disjoining pressures cause a maximum of Mc 
for d < dc but have nearly no influence for d > dc (dashed lines). 

Fig.[7| Maxwell constructions (horizontal lines) based on the local energy are given for (a) 
three different system thicknesses d (see legend) of a Rayleigh- Taylor unstable system 
at G = 5, C-i = 20, Hi = H2 = 0.01, M = 0, p = 1.826 and p = 0.1826. For 
d = 2, the Maxwell point is 1 (vertical dotted line), d > 2 leads to drop solutions 
(dashed line). Hole solutions are expected for d < 2 (solid hne). (b) Illustrates the 
occurrence of metastable states using G = 10, Hi = 0.1, H2 = 0.01, p = 1.826, and 
p = 0.1826, i.e. an Rayleigh- Taylor unstable system. The isothermal case (dashed 
line) is linearly unstable, whereas heating from below (M = 10, solid line) stabilizes 
the system linearly. However, the local maximum still exists at /i < 1 indicating a 
metastable flat film. 

Fig.lSl Shown is the critical Marangoni number Mc versus system thickness d for G = 10, 
C-i = 20, p = 1.826, p = 0.1826, A = 0.598, Hi = 0.1 and H2 = 0.01. The isothermal 
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system is Rayleigh- Taylor unstable (M = 0). For d > {d < dc) heating from above 
(below) stabilizes the system. However, for the shown range of M the system remains 
metastable, i.e. finite disturbances larger than a critical nucleus will grow. 

Fig. El Snapshots from of the nonlinear evolution of the interface for d = 1.3, G = 10, 
C-^ = 20, At = 0.1, Ax = 0.5, p = 1.826, /i = 0.1826, A = 0.598, Hi = 0.1 and 
H2 = 0.01. Without heating (M = 0) the system is linearly Rayleigh- Taylor unstable 
implying the growth of infinitely small disturbances (dashed line). For M > M^. = 4.62 
the system is metastable. We applied for M = 10 > Mc a strong disturbance of the 
interface (/;,«i0.7±0.1)in the vicinity of x = 50. 

Fig. El Three-dimensional snapshots of the long-time evolution of a Rayleigh- Taylor unstable 
system at d = 3, G = 5, G'^ = 20, p = 1.826, /i = 0.1826 and Hi = H2 = 0.01. The 
system size is = Ly = 200 with a resolution At = 0.1, Ax = Ay = 2. Initial (small) 
perturbations of the fiat interface lead to drop formation, and subsequent long-time 
coarsening. Finally, one single drop survives (not shown). 

Fig. ^2 Three-dimensional snapshots of the long-time evolution of a Rayleigh- Taylor unstable 
system at d = 1.3, G = 20, G'^ = 20, p = 1.826, p = 0.1826 and Hi = H2 = 0.01. 
The system size is L^. = = 100 with a resolution At = 0.1, Ax = Ay = 0.5. From 
initially small perturbations holes start to evolve [t = 150) smoothly. Subsequently 
long-time coarsening sets in at t ~ 1000 and finally an one-hole solution is reached 
(t > 2 X 10^). 

Fig. El Three-dimensional snapshots of the long-time evolution of a maze structure in a 
Rayleigh-Taylor unstable system at = 2, G = 5, G'^ = 20, p = 1.826, p = 0.1826 
and Hi = H2 = 0.01. The system size is = Ly = 200 with a resolution At = 0.1, 
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Ax = Ay = 2. It is clearly visible that neither holes nor drops are energetically 
preferred. In the long-time evolution {t > 1000) the usual coarsening takes places. 

Fig. El Shown are (a) the mean wavenumber (k) and (b) the energy E in dependence on time 
for Rayleigh- Taylor unstable systems with p = 1.812, fi = 0.1826, Hi = H2 = 0.01, 
and different thicknesses d as given in the legend. Horizontal thin lines give the 
corresponding fastest linear wave numbers. 

Fig-^S Grey-level plots of the interface height h at two timesteps and the difference of the 
two images for the time evolution in Fig.^J In the difference plot dark (light) areas 
indicate mass loss (gain) of the lower layer (a) During the first stage of the long- 
time evolution neighboring drops move towards each other to merge indicating the 
dominance of the translational mode of coarsening, (b) At a later stage, small drops 
shrink and neighboring large drops grow, indicating the dominance of the mass transfer 
mode of coarsening. 

Fig-El Given are (left) film thickness profiles h{x, t) and (right) x-components of the fluid 
velocity u{x = 20, z,t) for G = 5, C'^ = 20, p = 1.826, p = 0.1826, Hi = H2 = 
and a resolution At = 0.1, Ax = 0.5. The left panels also show contour lines of the 
streamfunction ip {u = —dzf)- (a) In a system where p < pc {d = 3) at t = 120 the 
interface height h{20,t) < 1 and dxh{20,t) < provides m < at the interface. The 
direction of the interface velocity corresponds to the solid arrows in Fig.|Sl(a), i.e. the 
lower layer is the driving layer, (b) If p > pc {d = 1.3), at t = 1200 the interface 
velocity u at x = 20 is positive and the upper layer is the driving layer (compare to 
solid arrows in Fig.EKb)). 

Fig. El Three-dimensional snapshots of the long-time evolution of a Marangoni instability for 
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M = -10, d = 2, G = 5, C-^ = 20, p = 0.548, /i = 5.48, A = 1.671, Hi = 0.01 
and H2 = 0.05. The system size is = Ly = 100 with a resolution At = 0.1, 
Ax = Ay = 1. From initially small perturbations the systems evolves smoothly to a 
single drop solution in the long-time limit (t > 10^, not shown). 

Fig-El Shown are (a) the mean wavenumber (k) and (b) the energy E in dependence on time 
for two thermocapillary and one electrohydrodynamic unstable system (see legend) 
with p = 0.548, /i = 5.48 and A = 1.671. Horizontal thin lines give the corresponding 
fastest linear wave numbers. 

Fig. El Three-dimensional snapshots of the long-time evolution of a Marangoni instability for 
M = 70, = 4, G = 5, C-i = 20, p = 0.548, p = 5.48, A = 1.671, Hi = 0.05 
and H2 = 0.01. The system size is = Ly = 200 with a resolution At = 0.1, 
Ax = Ay = 1. We started with initially small perturbations. At t ^ 700 one hole 
starts to evolve rapidly and subsequently more and more holes arise. At t ~ 1 100 long- 
time coarsening sets in and continues until a single large hole is reached (t > 2 x 10^). 

Fig-CHI Shown are the rescaled mobilities Q'l = G{l—p) Qi and Q2 = M Q2 for thermocapillary 
unstable systems with c? = 2, M = 70 (thin lines, numerical run in Fig.ll6|) and d = 
4, M = —10 (thick lines, numerical run in Fig.ll8|). For d = 4 the zero crossing of Q2 
is at h ^ 1.2. This leads to a suppressed interface evolution ior h > 1 resulting in a 
rapid hole evolution. For d = 2 smooth drop evolution results since the regions h > 1 
and h < 1 are roughly symmetric. 

Fig. EDI Three-dimensional snapshots of the long-time evolution of an electrohydrodynamic 
instability for d = 4, C'^ = 20, U = 30, e = 1.3, p = 5.48 and Hi = H2 = 0.01. The 
system size is = Ly = 300 with a resolution At = 0.1, Ax = Ay = 3. One finds a 
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smooth short-time evolution of drops. The long-time coarsening sets in at i 10000 
and the long-time scaling exponent is very small. 



48 



^0 




substrate 2 

/////// 



liquid 2 (^^2' P2) 



liquid 1 (l^pPi) 



h(x,y,t) 



\\\\\\\\\\\\ 

substrate 1 



FIG. 1: D.Merkt, Physics of Fluids 



49 




FIG. 2: D.Merkt, Physics of Fluids 
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FIG. 6: D.Merkt, Physics of Fluids 
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FIG. 7: D.Merkt, Physics of Fluids 
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FIG. 8: D.Merkt, Physics of Fluids 
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FIG. 10: D.Merkt, Physics of Fluids 
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FIG. 12: D.Merkt, Physics of Fluids 
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FIG. 13: D.Merkt, Physics of Fluids 
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FIG. 17: D.Merkt, Physics of Fluids 
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FIG. 18: D.Merkt, Physics of Fluids 
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FIG. 19: D.Merkt, Physics of Fluids 
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